Parametric Projection Operator Technique for Second Order Non-linear Response 
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We demonstrate the application of the recently introduced parametric projector operator tech- 
nique to a calculation of the second order non-linear optical response of a multilevel molecular 
system. We derive a parametric quantum master equation (QME) for the time evolution of the 
excited state of an excitonic system after excitation by the first two pulses in the usual spectro- 
scopic four-wave-mixing scheme. This master equation differs from the usual QME by a correction 
term which depends on the delay r between the pulses. In the presence of environmental degrees 
of freedom with finite bath correlation time and in the presence of intramolecular vibrations we 
find distinct dynamics of both the excite state populations and the electronic coherence for different 
delays r. 



I. INTRODUCTION 

Recent years have seen a rapid experimental develop- 
ment in multidimensional coherent spectroscopy. Devel- 
oped first in the nuclear magnetic resonance |l| it has 
been brought to infra red [2H4J and later into near infra- 
red and optical domains 0, Q . Since then it yielded new 
insights into photo-induced dynamics of electronic ex- 
cited states of small molecule s B , polymers [H, large 
photosynthetic aggregates [l, [Iol | and even solid state 
systems The two-dimensional (2D) Fourier trans- 

formed spectrum completely characterizes the third or- 
der non-linear response of a molecular ensemble in am- 
plitude and phase [T3 . Il3| providing thus, the maximal 
information accessible in a three pulse experiment. In 
the three pulse degenerate four wave mixing (FWM) ex- 
periment [14j, two independent adjustable parameters - 
the delays between the three pulses - allow experimental- 
ists to address the time and delay dependent third order 
non-linear response of the molecular system S^'it, T, r) 
via complete characterization of the third order non- 
linear signal Ei 3 \t,T,r) (see Fig. [I]). This informa- 
tion is conveniently represented by a 2D plot of the 
Fourier transformed signal E$ (wt,T, oj t ) which corre- 
lates the absorption frequency uj t with the frequency w t 
of stimulated emission, ground state bleaching and ex- 
cited state absorption contributions, separated from the 
absorption event by an adjustable waiting time T. Fem- 
tosecond photo-induced evolution of molecular systems is 
reflected in the amplitude, position and shape of 2D spec- 
tral features as functions of the waiting time revealing 
thus important information about the molecular struc- 
ture, intramolecular dynamics, as well as the interaction 
of molecular electronic transitions with their environment 

By analyzing 2D spectra of photosynthetic Fenna- 
Mathews- Olson (FMO) protein, Brixner et al. have first 
demonstrated the ability of the new technique to improve 
our knowledge of e nerg y transfer pathways in photosyn- 
thetic proteins |l5l . [18]. Later, signatures of electronic 
quantum coherence in this system have been theoreti- 
cally predicted [l6[ and experimentally verified |9]. The 



same measurements led to a surprising discovery of the 
long life time of these electronic coherence and to an iden- 
tification of a host of other effects which contribute to the 
wave-like nature of the energy transfer in these systems. 
Recently, coherent oscillations in 2D spectra of photosyn- 
thetic systems were demonstrated even for room temper- 
ature [l^ . [20| . The impact of these new finds has been 
felt well beyond the photosynthesis research community, 
and the quantum properties of energy transfer in bio- 
logical systems have attracted researches from seemingly 
unrelated fields of quantum computation and quantum 
information science. Questions about the relevance of 
quantum effects in natural light harvestin g ha ve led re- 
searches to study quantum entanglement [2ll423| . vari- 
ous aspects of the environmental assistance in quantum 
transport and optimality of transport processes (23-[28l|. 
Meanwhile, measurement of electronic coherence have 
found utility in determining structure-related properties 
of photosynthetic systems. Thus, the electronic beating 
of 2D spectra have been recently used to precisely de- 
termine electronic energy levels in congested spectra of 
molecular aggregates [29j . 

In order to yield the above discussed impressive re- 
sults, optical 2D experiments have to be accompanied 
by a thorough analysis, which requires detailed theoret- 
ical understanding of the molecular response to excit- 
ing light. For the description of the most 2D experi- 
ments, the third order semi-classical light-matter inter- 
action response function theory is well established. Re- 
sponse functions of model few-level systems with pure 
dephasing (i.e. with no energy transfer between the lev- 
els) can even be expressed analytically in terms of the so- 
called energy gap correlation function (EGCF), using the 
second order cumulant in Magnus expansion |14j . Some 
examples of small chromophores in solution fall in this 
category [3(ill3l| when investigated on time scales shorter 
then radiative life time. For a Gaussian bath this ana- 
lytical theory is exact, and thus knowing the EGCF of 
the electronic transitions enables us to determine linear 
(absorption) as well as non-linear spectra. 

However, the construction of exact response func- 
tions for realistic energy transferring systems, such as 
Frenkel excitons in photosynthetic aggregates [HI, is 
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no longer possible. Photosynthetic complexes are rela- 
tively large, and the proper methods to simulate finite 
timescale stochastic fluctuations at finite temperatures 
[UHlfl carry a substantial numerical cost. Practical 
calculations thus require some type of reduced dynam- 
ics where only electronic degrees of freedom (DOF) are 
treated explicitly. These approaches usually rely on a 
host of approximations that seem to work well for most 
spectroscopic techniques [36 39] and in systems where 
the details of system-bath coupling are apparently less 
important With increasing details of the excited 

state dynamics revealed by the 2D spectroscopy [4ll - |43j 
and with increasing size of the studied systems, it be- 
comes more and more important to keep the numerical 
cost of simulations low, while simultaneously account for 
experimentally observed quantum effects. One of the 
possible research directions is to extend on existing re- 
duced density matrix (RDM) theories Q EH or relax 
certain approximations [46|. 

It was demonstrated that RDM master equations de- 
rived by projection operator technique reproduce exactly 
the linear response [47]. In the case of higher order re- 
sponse functions, however, the same approach necessarily 
neglects bath correlations between the different periods 
of photo-induced system evolution (i.e. between so-called 
coherence time r and the waiting time T). The fail- 
ure to account for this correlation leads sometimes to a 
complete loss of the experimentally observed dynamics 
in simulated 2D spectra, such as in the case of the vibra- 
tional modulation of electronic 2D spectra [30(. Because 
vibrational modulation leads to effects similar to those 
attributed to electronic coherence, developing methods 
that can account for its effect reliably in complex sys- 
tems is of utmost importance. One possible approach 
to the problem is to derive equations of motion for the 
response as a whole, and to take the previous time evo- 
lution of the system into account explicitly (4§j . In this 
paper, we take a different route in which the correlation 
effects are treated by a specific choice of the projection 
operator. Once a projector is specified, it can be used 
for any method of treating system bath interaction. We 
apply the previously suggested parametric projection op- 
erator technique [49| to a calculation of the second order 
response of a quantum system. The second order re- 
sponse operators can be used to determine the state of 
a molecular system subject to excitation by a weak light 
with arbitrary properties, and thus its importance goes 
beyond the semi-classical system light interaction the- 
ory [5(| • Alternatively, calculation of the second order 
response can be viewed as a first step towards more in- 
volved calculation of the third order response functions 
that is required for the modeling of the third order non- 
linear spectra. 

The paper is organized as follows: in the next section 
we specify the model Hamiltonian and the description 
of the system's interaction with the light and its envi- 
ronment (the thermodynamic bath) . In Section IHII we 
discuss the third order response functions of a multilevel 



systems, and we point out general limitations of their 
evaluation by the reduced density matrix master equa- 
tions. In Section lPvl we write out the density operator de- 
scribing the excited state of a molecular system in terms 
of the second order response operator, and we discuss the 
advantages of usage of the parametric projection opera- 
tor over the standard projectors. The details of the para- 
metric projection operator and the corresponding master 
equation are introduced in Section [V] Numerical results 
for the excited state dynamics are then discussed in Sec- 
tion [VTJ Some definitions and calculation details can be 
found in the Appendices. 



II. MODEL SYSTEM AND BLOCK 
FORMALISM 

Our primary interest lies in photosynthetic aggregates 
of chlorophylls. Frenkel exciton model is well established 
for the description of photosynthetic systems and their 
spectroscopy [i~5l |32||. The second order response, which 
we have in mind here, corresponds to the interaction of 
the excitonic system with first two pulses in a standard 
FWM spectroscopy (see Fig. [T]), and it should corre- 
spondingly describe the time evolution of the system in 
the excited or ground state. The part of the response 
associated with evolution in the excited state also corre- 
sponds to the material quantities that govern excitation 
of molecular systems by arbitrary quantum light [50] . 
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Figure 1: A scheme of a typical FWM experiment. Three laser 
pulses with delays r and T are incident on a sample. A non- 
linear signal generated into a selected direction, — fci +k2 + k$ 
here, is detected. 



For a second order response we need to consider only 
single exciton band, i.e. the collective excited states of 
the aggregate with only a single aggregate member ex- 
cited. The Frenkel Hamiltonian including environmental 
contributions has the general form H — Hs+Hb+Hs-b, 
where the purely electronic (system) Hamiltonian H$ 
reads 

H S = e» \g)(g\ + ]T (e n + (K - &)) \n)(n\ 
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the bath Hamiltonian Hb can be identified with 



and the system-bath interaction H,s-b reads as 
H s - B = ^2^n\n)(n\. 



(1) 



(2) 



(3) 



Here, \g) is the electronic ground state of the aggregate of 
molecules with energy e 9 , \m) are electronic states with 
the m th molecule of the aggregate of N molecules excited 
with energy e n , j nm represents resonance coupling be- 
tween excited states \m) and |n), <& represents potential 
energy surfaces (PES) of the bath DOF in corresponding 
electronic states, and T is the kinetic energy of the bath 
DOF. We also defined so-called energy gap operator 



A3„ = 3« -#»-(#«- 



(4) 



where (•) represents the equilibrium quantum mechan- 
ical average over the bath DOF. Symbol • denotes an 
arbitrary operator. The interaction of the system with 
light will be described by the usual semi-classical light- 
matter interaction term in dipole approximation 



H s - E (t) = -M-E(i). 



(5) 



Here, /j, is the transition dipole moment operator of the 
aggregate system and E(t) = eE(t) is the electric field 
vector of the light with polarization vector e. 

The states \m) with excitations localized on individual 
member molecules of the aggregate allow us to define the 
properties of the aggregate based on crystal structure 
and quantum chemical inputs. It is, however, often more 
convenient to work in a basis of the eigenstate |m) of the 
electronic Hamiltonian Hg. For later convenience we will 
define projection operators K m — \m){m\ projecting on 
a state in the local basis and K m — \rh){rh\ projecting on 
the eigenstate (excitonic) basis. For the discussion of the 
non-linear response functions, it is advantageous to use a 
superoperator formalism. We define the dipole moment 
superoperator as 



and the Liouvillian as 



£• = -[#,.]_. 

We will also use the evolution operator 



E/(*)=exp 

and the evolution superoperator 

U{t)» = U(t)»U*(t), 



(G) 



(7) 



(8) 



(9) 



where convenient. 

In the second order system-bath coupling theories, the 
energy gap operators, Eq. (j4}, are present in the equa- 
tions via the EGCFs 



C mn (t) = (Ui(t)A$ m U(t)A$ n ). 



(10) 



Double integration of C m n(t) over time yields so-called 
lineshape functions 



g mn (t) = / dl ~ I dl ~' C« m ( T ')> (ID 



which determine absorption and emission line shapes. 
Later in the text, it will be convenient to define the en- 
ergy gap operators in the excitonic basis 



A$ ai = J2^k(&\k){k\b), 



the excitonic EGCFs 



C- a i(t) = (tf(t)A$ m U(t)A<Z> u ) 

= Y / \(m\ 2 \(W)\ 2 c kl (t), 



(12) 



(13) 
(14) 



kl 



and the excitonic lineshape functions 



9- a l(t) 



ft 2 



dr / dr' C ai (r'). 



(15) 



We are interested in the influence of the bath DOF 
on the dynamics of the electronic DOF in the excited 
state. The bath can consists of several types of collec- 
tive modes which can be described by various types of 
bath correlation function. Experimental situation is typi- 
cally well described by one or several overdamped Brown- 
ian oscillator modes standing for an macroscopic number 
of harmonic DOF, and several underdamped oscillator 
modes standing for some important vibrational coordi- 
nates, e.g. normal modes of the chromophore molecules 
[30j . Both of these limits can be conveniently described 
by by the gen eral Brownian oscillator correlation function 
C{t) (see [1J|) which is a function of temperature T, fre- 
quency f^Bathof the oscillator, damping coefficient 7 and 
the reorganization energy A. In the case of a strongly 
overdamped mode i.e. when 7 3> 2ilbath the formula 
simplifies significantly. At high temperatures we have 

C(t) = fi.AAcot(A/i/3/2)exp(-A<) 



ihXA exp(-Ai), 



(16) 



with A = I/7 . In the opposite case of a non-damped 
oscillator 7-^0 the correlation function reads 

C(t) = Af2 Ba thft COth(/3fiSl B ath/2)cOs(riBathi) 
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left to right) follow the double-sided Feynman diagram 



i sin(0 Ba th*) 



(17) 



The Hamiltonian, Eq. (p}, has a block form. In case 
we would consider excited states with up to N excita- 
tions per aggregate a more general Hamiltonian would 
have to be written including N + 1 blocks with one block 
containing just the ground state \g), and the N blocks 
containing the excited states with one, two, . . . , up to 
N excitations. The block structure is enabled by the 
fact that the Hamiltonian operators, Eqs. fT]) to ([3j). 
do not contain any terms that enable de-excitation of 
excited states. This is a well justified assumption in 
studies of ultra-fast dynamics of the chlorophyll based 
photosynthetic aggregates [5l|. The transitions between 
different blocks of the Hamiltonian are only enabled by 
the interaction with the light, Eq. ([5]). The second or- 
der response, which is of interest in this paper, requires 
on the ground- and the single-exciton blocks. We will 
therefore use upper indices g (ground state) and e (sin- 
gle exciton block) to denote different blocks of operators 
and superoperators whenever we are interested in opera- 
tions on and evolution of their individual blocks. Thus, 
e.g. the time evolution of the system in excited state 
single exciton manifold is described by the RDM block 
P ( ee \t) = U^ eeee \t) P ^{Q) = U (ee \t)^ ee \0)U^(t), 
and the action of the dipole operator (superoperator) pro- 
motes the ground state block into a coherence block as 
p( e 9)(t) = Hapten) (t) = V {e999) p {99 \t). When working 
in a particular basis of states (e.g. {|m}} = {|1), |2), . . . }) 
we can write out the sums over the states explicitly, e.g. 
as 



i run 



w = E^ } w^ e) (o) 



kl 



(18) 



A7 



The rules are rooted in the simple well-known fact that 
matrices can be multiplied by blocks. When discussing 
the most common non-linear experiments, the block 
structure has to be considered up to the two-exciton 
block. 



III. MULTI-POINT CORRELATION 
FUNCTIONS IN NON-LINEAR SPECTROSCOPY 

Non-linear response functions have in general the form 
of multi-point correlation functions. For a two band sys- 
tem (ground state and single excitons), the third order 
non-linear spectroscopy is completely described by four 
response functions [14], listed in Appendix |Aj As an ex- 
ample, we will consider the response usually denoted as 
i?2, Eq. (|A2[) . We can notice that the block (upper) in- 
dices of the evolution operators in Eq. (|A2[) (read from 



in Fig. (see [14j for details). During the so-called 



population interval T of the response, the system evolves 
in the excited state band (ee) . If no resonance coupling 
is present {j mn — in Eq. ([!"])) each response function 
can be split into a sum of independent components. For 
i?2 this means 



R 2 (t,T,T)=J2 R 2,ki(t^ T l 



(19) 



ki 



where the components 

x *b$& (T)UiZ 9e) (r)W^ } (20) 



can be evaluated analytically in terms of the line shape 



functions, Eq. (fTl) . Here, W ( 



(gg) 

oq 



>cq\g)(g\ is the to- 
tal equilibrium density operator (the electronic energy 
gap is assumed to be much higher than k B T) and w eq 
represents the equilibrium density operator of the bath 
alone. For j mn > no analytical result is available, and 
one needs to resort to some master equation simulations 
of the evolution superoperators IA. However, the master 
equations can only deliver certain reduced (i.e. averaged 
over the state of the bath) version of this superoperator. 

The usual way of deriving master equations for the 
RDM is to apply a projection operator V which reduces 
the full density matrix W(t) of the system to its selected 
part [52], in our case, to the electronic DOF 



p(t) =tr B {W(t)}. 



(21) 



The most popular prescription, the so-called Argyres- 
Kelly (AK) projector (53|, reads as 



V W(t) = p(t)w, 



eq 



It is easy to see that inserting the identity 1 = V 
where Q=l-V into Eq. (HJ) leads to 



(22) 
+ Q, 

p(t) = ti B {U(t)VW(0) + U(t)QW(0)}. (23) 

Provided that QW(0) — (a condition satisfied in the 
first interval of the response function, Eq. (|A2[) ) one 
can define the reduced evolution superoperator U(t) = 
trB{U{t)w cq } such that p(t) — U(t)p(0). The evolution 
superoperator U(t) can be calculated from a master equa- 
tion that can be derived by projection operator formalism 

Let us apply the same method to higher order response 
functions. Let use assume we have derived a master equa- 
tion by applying the corresponding projector operator V 
and calculated elements of the reduced evolution super- 
operator U(t). We can then assemble an approximate 
response function 
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'^m\T)ufj e \r)p^{Q). (24) 



Here, we set all transition dipole moment elements to one 
for the sake of brevity. It can be shown that R 2t u can 
also be written as 



R 



xruir k r\T)r^ g r g r>(r) P 

However, the exact expression for R 2t ki reads as 
R 2 , k l(t,T,r) = tv B {uiHf(t)(V + Q)uiZr\T) 



,{gege) (n _\ n {gg) 



(0)}. 



(25) 



x (V 



(26) 



Eqs. (|25|) and (f26f differ by Q-containing terms that 
cannot be in general eliminated, and consequently one 
cannot expect master equations based on a single pro- 
jector operator V of any type to reproduce the third or- 
der response functions. This applies also to the validity 
of non-perturbative schemes of calculations of non-linear 
response, such as those derived in Refs. [54[ and [55) . 

A general solution of this problem was proposed in 
Ref. [49]. It was argued that one cannot write down 
a single exact master equation for all three intervals of 
the response function. Rather, one has to write a differ- 
ent master equation for each interval. This is formally 
possible by introducing three different projectors Vo (i-e. 
standard AK projector) for the first, P T for the second 
and P t +t for the third interval of the response [4^|. The 
projectors P T and P t +t are constructed so as to cancel 
the Q-containing term exactly for j mn — 0. In this limit, 
all response functions, Eqs. (|A1[) to (IA4[) . can be exactly 
reproduced by the corresponding master equations. In 
the following sections, we will treat the case of j mn ^ 0, 
for the case of the second order response, i.e. we will de- 
rive equations of motion for the reduced density matrix 
using the projector P T . The application of this approach 
to the full response function will be treated elsewhere. 
An alternative to the projector approach is to attempt 
to derive specific equations of motion for each response 
function as a whole in a specific perturbation scheme, 
e.g. second order convolutionless QME (CL-QME) as in 
Ref. [481] . While both approaches should lead to similar 
results, the projection operator technique is not limited 
to any specific way of expanding the equations of motion 
in terms of system-bath coupling, and we believe it is 
therefore somewhat more flexible. 



IV. SECOND ORDER RESPONSE TO LIGHT 

The second order optical response can yield an optical 
signal on the sum or difference frequencies of the exciting 



field. Here, we consider only the degenerate experiment 
when the excitation field have the same frequency, and 
we are interested in particular in the zero frequency part 
of the response which does not generate an optical sig- 
nal. Unlike in spectroscopy where we study a quantum 
mechanical expectation value (polarization or field), here 
we would like to study the state of the system achieved 
by the excitation. Correspondingly, the second order re- 
sponse will not be expressed in terms of response func- 
tion, but rather in terms of some response operators. By 
the "system" we mean electronic DOF, and thus the re- 
sponse operators correspond to the reduced density ma- 
trix of the system in the same way the response function 
corresponds to an expectation value of e.g. polarization. 



A. Excited State of a Molecular System 

Assuming that we excited the system by some external 
field E(t), the second order reduced density operator of 
the system reads as 

oo oo 

P {2 \t) = - f dti f dt 2 tT B {U(t 2 )V 



x U{tx)VW^\g){g\}E{t ~ t 2 )E(t - t 2 - h). (27) 

The excited state part of the second order density oper- 
ator p( 2 ^ then yields 



OO OO 

dti J dt 2 
o o 



Ri(t 2 ,ti)A*(t - t 2 )A(t - t 2 - ti)e 



%UJt\ 



+ R n {t 2 ,t{)A{t - t 2 )A*(t -t 2 - h)e- l ^\ , (28) 
where 

R I (t 27 t 1 )=Tr B {u^(t 2 ) 
x v (eee9) ^ (e9e9) (ti)V (e " 9) «;cq|<7)(.g|}, (29) 



and 



R 



ii 



(t 2 ,t 1 ) = Tr B {^ eeee \t 2 ) 



x V (eese) W (9ese) (ii)V (9e99) w cq |<7)(. 9 |}, (30) 



are the second order response operators. In Eq. (1281) we 
introduced the electric field envelopes A(t) and the car- 
rier frequency w so that E(t) — A(t)e~ lut + c.c. The 
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interaction of the molecular system with arbitrary light 
can be expressed using the operators Rj and Rji as dis- 
cussed in Ref. [5(|. For a general state of the light \ip ia d), 
the terms A* (t - 1 2 ) A(t - 1 2 -ti)e iujtl have to be replaced 
by the light correlation function I(t — t 2 ,t — t 2 — ti) = 
(Aad\E^(t-t 2 )E(t-t 2 -ti)|^ r ad)i where E(t) is the op- 
erator of the electric field in Heisenberg representation. 
This allows us to calculate the state of a system excited 
by light with arbitrary properties. 

In this paper, we will concentrate on the situation when 
the molecular system is excited by two ultra-short laser 
pulses traveling in two different directions k\ and k 2 of 
which the second one arrives with a delay r. Thus we 
have Ai(t) = A e tk i- r 5(t + T), A 2 (t) = A e ik ^ r 5{t). The 
time zero is set to the center of the second pulse. This 
corresponds to a typical situation in the coherent non- 
linear spectroscopy, and the time t in which the excited 
state evolves corresponds to the population time T of 
the non-linear spectroscopy. Inserting the delta function 
envelopes into Eq. ([2"5f and assuming that the pulse with 
wave- vector fci precedes the pulse with wave- vector k 2 , 
we arrive at 

pf\t;r) = R I (t,r)\A \ 2 . (31) 

This operator corresponds to the excited state time evo- 
lution in the left hand side (l.h.s.) diagram of Fig. [5J3. 
The base of this diagram corresponds to the so-called 
rephasing pathways. When the pulse sequence is reverted 
(pulse fci arrives second with delay r) we obtain so-called 
non-rephasing pathways and the time evolution corre- 
sponds to the operator 

p^{t;r)=Ru(t,T)\A \ 2 . (32) 

It is easy to verify that pj(t;r) and pf (i;r) do not 
have to be Hermitian, and they alone cannot be said to 
represent a state of a molecular system. This is the result 
of them being only a portion of the perturbation expan- 
sion of the non-linear response operator. Only the sum 
of Eqs. ([2l?j) and (|3"0|) yields a Hermitian operator which 
can describe excited state of a molecular system. This 
situation has an analogy in the case of 2D coherent spec- 
troscopy where the rephasing and non-rephasing signals 
alone cannot be interpreted as representing absorption 
and emission events, while the sum spectrum can be as- 
signed this interpretation [13j . When the system is ex- 
cited by a single finite length pulse, the two contributions 
to the excited state are equally weighted, guaranteeing 
thus the proper properties of the corresponding density 
matrix. 

In the following section we will discuss the equations 

(2) (2) 

of motion for p\ (t;r) and pjj (t; t) to asses the influ- 
ence of the delay r on the dynamics of their diagonal 
(populations) and off-diagonal (coherence) elements. We 
will thus attempt to answer the question whether in the 
non-linear spectroscopic methods, such as 2D coherent 
spectroscopy, we observe the expected excited state dy- 
namics, i.e. the one unaffected by the delay r. 



B. Simulation of the Response by Reduced Density 
Matrix Equations 

In order to calculate the second order non-linear re- 
sponse, we need to evaluate the component expressions 
of Eq. (|30p . We start with equations of motion for the 
perturbation of the total density matrix W(t). We define 
the first order operator as 

W$(t) = U^ e Ht)V^w cq \g)(g\. (33) 

This operator corresponds to the evolution after the first 
interaction of the system with light in response function, 
Eq. ([31?)) . The response function has to be read from 
left to right. The operator Wfj(t) satisfies the following 
equation 

^w£\t) = -i£^W$(t), (34) 

with the initial condition Wj/ 1 (t — 0) = 
V^ e ^w cq \g)(g\ = w eq \g){g\^ e \ RDM equation 
of motion can be found applying standard AK projec- 
tion operator, Eq. ([22]), for which (1 - Vo)W$(t = 0) is 
equal to zero identically. The procedure of the derivation 
of the RDM equation of motion leads to 

= ^r^pfiit) - n { r e) it) P fht), 05) 

where C^ 9 ^ is the coherence block of the electronic 
Liouville superoperator, and TZg 9ege ^ (£) is some second 
order dephasing tensor. The particular form of TZ^ 9e9e ^ 
depends on the approximations and the theory applied. 
For a pure dephasing model and a harmonic bath, the 
dephasing tensor which follows from a second order CL- 
QME (see e.g. [46]) can be shown to yield an exact dy- 
namics for pf] it) (4?| . In secular approximation, Eq. 
([55) yields 

(36) 

which is a set of independent equations for optical coher- 
ences. Note that we use index k to denote the excitonic 
representation. In the limit of j mn — > (pure dephasing) 
the states \k) and coincide. 

The evolution in the second propagation interval can 
be expressed through the density operator 

Wfftt-r) =U {eeee) (t)V ieeeg) W i I 1 I ) (T), (37) 
which satisfies 

§-Wlf(t; r) = -tC^W^it; r), (38) 

(2) 

with the initial condition Wj/ (t = 0; r) = 
yOese)^!) (yj _ fj,(eg) W (^ T y Thc app li ca tion of the 
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AK projector would lead to a loss of information in 
the construction of the response function, because the 
term Qop^W^ (t) = (1 - Vo)W^{r) ^ 0. Con- 
sequently, even if we successfully derive an exact mas- 

(21 (2) 

ter equation for p\j{t) — Ttb{W u (t)}, the response 
function constructed from this equation would be miss- 
ing the Q-containing terms. The parametric projec- 
tion operator recently suggested in Ref. [49| has the 
property of eliminating the initial term approximately, 

(1 — "Pt) fJ-egWii '(r) ~ . In the case of pure dephasing, 
it turns to zero exactly. We can therefore write an equa- 
tion of motion for pfj which is analogical to Eq. (|35|) 



n^ ee) {t)pf){t;T). 



(39) 



Again, the £^ eeee ^ is the electronic Liouville superoper- 
ator and T?jf eee ^ (t) is the superoperator describing the 
electronic energy relaxation and the electronic coherence 
dephasing in the single-exciton band. In the following 
section, we will derive the relaxation tensor f^. eeee ^ (£j 
corresponding to the projection operator P T . 



V. MASTER EQUATIONS WITH 
PARAMETRIC PROJECTOR 

In this section, we will apply the well-known Nakajima- 
Zwanzig identity in the second order in system bath in- 
teraction Hamiltonian together with the parametric pro- 
jection operator proposed in Ref. [49]]. We are interested 
in the second interval of the non-linear response, and in 
the evolution in the single exciton band in particular . 



A. The Choice of Projector 

The projector V r is chosen in such way that it contains 
time evolution of the bath in the first interval of the re- 
sponse, where the relevant system dynamics is the one of 
an optical coherence. The projector corresponding to the 
pathway Ru and the length of the first interval r reads 
according to Ref. [49] 



V T ,II* 



tr B {«} C/ 9 (r)w eq J2 U n(r)K n e g '^ . (40) 



This choice gives an exact description of the bath for 
jmn = 0. For a non-zero coupling it corresponds to the 
secular approximation in the first interval equation of 
motion, Eq. (|3"6"]h Further on in the text, all derivations 
will be done for the pathway Ru, and we omit the index 
77. The treatment of the pathway Ri is analogical. We 
define evolution operators describing the dynamics of the 
environmental DOF while the system is in its electronic 
ground state 



C/ 9 (r)=exp(--(T + $ 9 )r), 



and in the excited eigenstate \h) 



(41) 



C/?(r) = exp l--(T + *5 fi - - &))t \ . (42) 

By the choice of the projector, Eq. (|4U1) . we prescribe an 
ansatz 

W^(0) = W<»W(t) = p««(r)™«, (43) 



where 



(44) 



For zero resonance coupling jmn, this is an exact pre- 
scription for the bath. For non-zero resonance coupling, 
it is an approximation comparable to the secular approx- 
imation. Projector, Eq. (|4T)|) , can be written in short 
as 



V T * = tr B {.}i4 J 



It is important to note that wi is not a purely bath op- 
erator, and it does not generally commute with p. Also 
the interaction picture with respect to the system Hamil- 
tonian Hs denoted by (7) applies to it. It stands on the 
right hand side (r.h.s.) of p when evaluating Ru, while 
in Rj it would stand on the l.h.s. of p. This follows from 
the diagrams in Fig. The full form of W^(t) is 

W {1) {t) = U 9 (t)w c(1 x 

^^WJfBe^+W-^JWV^O). (46) 

h 

From now on, the upper index (2) will be omitted in text 
for the sake of brevity. 
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Figure 2: Double-sided Feynman diagrams of the third and 
the second order response. Part A: The Feynman diagram of 
the Liouville pathway 7?2. Part B: In full lines, the Feynman 
diagrams of the response operators Ru (left) and Ri (right). 
The dashed part completes the diagram into a corresponding 
third order response. 
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We can verify that the projector property V 2 = V T is 
fulfilled 



VrVr (\m)(h\ w) 

= \m)(n\ tr fl {U\r)w cq C#(t)} ti B {w} e 2 ^) 

= \m)(n\e- g '™ (T hr B {w} e 2ff ^ (T) = V T (\m)(n\ w) . 

(47) 

Here, we used expression for the line shape function g(t) 
in the second cumulant approximation 



e-^) =tv B {U^(t)U e (t)w eq } 



(48) 



The action of the projector V T (= V r ,ii) on the electronic 
state is asymmetric, because the projector was derived 
for the Liouville pathway II (see Fig 03). 



B. Parametric Master Equation 

Now, we apply the projector V T , Eq. (|45p. to the 
Nakajima-Zwanzig identity in the interaction picture 



—V T W {I) {t) = NZi + NZ 2 + NZ 3 = 

-V T C (I \t) exp (^- l J dT ' 2r/: (7) (r')Qrj QrW(t ) 

t~to r /i \ 

J dr' V T C (I \t) exp (-t J\t" Q t C {i) {t")Q t 



Q T C (I \t - T ')V T W {I \t - t') 

iVrC {I) (t)VrW {I) (t). 



(49) 



We will use Eq. (|49p up to the second order in D- 1 ', 
and we set to = 0. The term NZi corresponds to so- 
called initial term, which we made equal to zero by the 
choice of the projector. The last term NZ3 corresponds 
to an effective Liouvillian, and it is usually a purely elec- 
tronic operator. Now, with the parametric projector it 
contains additional terms originating from the system 
bath interaction. Its purely electronic part will stand 
for the effective Liouvillian jC,^ eee ^ in Eq. (|39f . while 



its additional r-depending contribution we will add to 
the relaxation tensor 7^t ee ' . Finally, the term NZ2 is 
a starting point for the derivation of the second order 
relaxation term, which has a form of a convolution be- 
tween the RDM and some memory kernel. In the second 
order expansion of Eq. (|49l) and with an approxima- 
tion W^(t — t') « W( T >(t), the resulting equation for 
V T W^ (t) coincides with the second order approximation 
of the equivalent time-convolutionless identity [56] . Thus 

the NZ2 will lead to a contribution to the tensor TZi eeee ^ 



in Eq. (|39p . However, the parametric projection opera- 
tor technique allows us to keep the convolution form of 
the equation if it is desired. While the convolution form 
may have some advantages over the CL-QME [46], only 
the CL-QME leads in the limit of j mn = to a result co- 
inciding with the one obtained by the second cumulant 
treatment of the non-linear response functions (see e.g. 
Rcf. HH and the Appendix [D} 



Let us point out the most important aspect of the 
application of the projector V T with the identity, Eq. 
(|4"9")l . It is important to note that the projector T T itself 
contains the system-bath interaction to all orders in the 
form of the exponential of the line-shape function (see 
Eq. (|4"tJ10 . The success of the second order master equa- 
tions (of the form p = —a^p, where the dot denotes the 
time derivative, and ai is some second order operator) 
lies in the fact that their solutions includes all orders 
of the perturbation (p — e~ a2t po). The solution corre- 
sponds to a partial summation of the perturbative series 
to infinity. For some types of bath, such as the bath 
consisting of harmonic oscillators, this may even lead to 
exact master equations [47]. When higher order terms 
are added to the right hand side of the equation motion 
by a procedure that does not respect the form of higher 
order terms dictated by the Nakajima-Zwanzig identity, 
the resulting equation of motion may lead to unphysical 
results. Therefore, one has to take care in application 
of the projector V T , not to allow higher than second or- 
der contributions to appear on the right hand side of Eq. 
(l4"9"|) . Since the difference between projectors Vo and V T 
is only in dynamics of the system-bath coupling during 
time r, their difference is at least of the first order in A$. 
Since the NZ2 with AK projector is already of the sec- 
ond order in A $ in all its terms, the difference caused by 
using projector V r will be of higher order in A$. In the 
second order master equation, the term NZ2 with pro- 
jector V T has to be equivalent to the form obtained from 
with AK projector (see e.g. (52j). Applying the following 
approximation p^\t — t') w P {t) m the term NZ2 we 
obtain it in the form the second order relaxation term of 
CL-QME 



The details of the evaluation of the term NZ3 are 
presented in Appendix [Bj As expected it yields a r- 
dependent term. Putting all results together and tracing 
over bath DOF we obtain 
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t 

mn q 

- 9mn{s)K m {t)K n (t - S)p {1 \t) 

+ g: nn (s)K m (t)p {I) (t)K n (t-s) 
(s)K n (t ;- s)pV\t)K m {t) 
-g; im ( S ) P (I) (t)K n (t- s )K m (t) 

+ J2 (K m (t)pW{t)K- k pV\t)K~ k K m {t)) 
m k 

x(g* mk (t + T)~g* mk (t)) . (50) 

The first five lines of Eq. ([50j) corresponds to the stan- 
dard CL-QME, while the last two lines represent the 
r-dependent contribution which the standard CL-QME 
does not predict. 

Let us investigate the r-dependent term only. First, 
we turn to Schrodinger picture by substituting p^ (t) = 
Ul(t)p(t)U s (t). We denote the new r-dependent term of 
the CL-QME by V(t; r) 

V(t; T )p(t) = 

£ (K m (t)p(t)K k -p(t)K k K m (t)) 

mk 

It can be easily verified that the new 
term preserves the trace of p(t), because 
Tr (K m {t)p{t)K k ~ p(t)K k K m {t)) = 0. ' From the 
inspection of the term g* ~{t + r) — g* ~(t) we can 
conclude that the effect of the parameter r is transient. 
If EGCF tends to zero on a time scale given by some 
bath correlation time, this term also tends to zero. 
The dynamics at long times is therefore not affected by 
the delay between the two excitation interactions. In 
Appendix [C] we derived the difference term, Eq. ([ST]) , 
of a special case of molecular homodimcr with bath 
fluctuations uncorrelated between the sides. We found 
the r-dependent term to be identically equal to zero in 
this case. In the next section we study the effect of the 
r-dependent term numerically for a heterodimer. 

VI. NUMERICAL RESULTS AND DISCUSSION 

In this section, we study the dynamics of the elements 
of the response functions, Eqs. $Z9§ and ([30]). The RDM 
p(t; r) for which we derived Eq. (|50|) corresponds the 
response function Rjj (see Eq. (|32[) for the case of A$ = 
1). We compare the dynamics in two cases. In the first 
case, the evolution operator £/( eeee ) (T) appearing in Eqs. 
(p?9")) and (|3"0"|) is calculated by standard time dependent 
CL-QME derived using AK projector, Eq. In the 

second case, it is calculated using Eq. (I5T)1) . We use 



the relation pi,ij(t) = p* u ^(t) between the RDMs from 
different Liouville pathways. 

The most simple system which exhibits r-dependent 
correction to the standard CL-QME is a molecular het- 
erodimer. In general, it is characterized by orientation 
and magnitude of its transition dipole moments, the ex- 
cited state energies e\, £2 of the component molecules, 
their resonance coupling J and the properties of the bath. 
We denote the difference of the excited state energies by 
A = £1 — £2, and we set the magnitudes of the transi- 
tion dipole moments to unity. The initial condition is 
assumed in a form p^(0) — \g)(g\. The evolution during 
the first interval of the response follows Eq. ([55)) . For 
the calculations presented on Figs. [31 [Hand [SI we choose 
the anti-parallel orientation of the transition dipole mo- 
ments, while for the calculation shown on Figs. [B]and[71 
we choose the parallel orientation. The temperature is 
set to T = 300 K in all calculations. 

First, let us investigate the sensitivity of the excited 




T[fs] 



Figure 3: The time evolution of the real parts of Ri and Ru 
for a heterodimer with energy gap A = 100 cm -1 , interacting 
with a bath represented by a general Brownian oscillator with 
parameters A = 50 cm -1 , 7 = 1 ps _1 and f^Bath = 100 ps -1 . 
Red lines, the matrix element 12 of Ri according to stan- 
dard CL-QME (dashed) and the corresponding parametric 
CL-QME (full). Black lines, the matrix element 12 of Ru ac- 
cording to standard CL-QME (dashed) and the corresponding 
parametric CL-QME (full). 
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Figure 4: The time evolution of the real parts of Ri and Rn 
for a heterodimer with energy gap A = 100 cm -1 , interacting 
with a bath represented by a general Brownian oscillator with 
parametersA = 5 cm" 1 , 7 = 1 ps" 1 and ^Bath = 19 ps -1 . Red 
lines, the matrix element 12 of Ri according to standard CL- 
QME (dashed) and the corresponding parametric CL-QME 
(full). Black lines, the matrix element 12 of Rn according to 
standard CL-QME (dashed) and the corresponding paramet- 
ric CL-QME (full). 



state dynamics to the interplay of the delay r and the 
phase of the bath vibrations. Fig. [3] shows the dy- 
namics of the electronic coherence between the excited 
states |ei) and e2 } in presence of the bath represented 
by a single-mode general Brownian oscillator with pa- 
rameters A = 50 cm" 1 , 7 = 1 ps" 1 , OBath = 100 ps" 1 and 
A = 100 cm" . Both calculations show that the stan- 
dard CL-QME calculation of ^( eeee ) (T) is insensitive to 
the phase of bath vibration mode during the time evolu- 
tion in the first interval. The parametric CL-QME, Eq. 



(1ST))) , shows a distinct sensitivity to this phase. In both 
cases, the resonance coupling J is chosen to be zero, and 
the dynamics can therefore be evaluated exactly by the 
cumulant expansion technique (see Appendix [Dj) . The 
numerical evaluation using Eq. (|50[) indeed matches the 
analytical result. In the calculation on the Fig. [31 the 
vibration of the bath is much faster than the period (333 
fs) of the electronic coherence. The two theories give 
the same result for r = fs, then they start to devi- 
ate and after one period of bath oscillator, at r = 60 fs, 
they coincide again. Initial phase of pi t \2 and pn^vi is 
in general different at T = fs because of their different 
time evolution in the first interval (they are not simply 
complex conjugates of each other). In the Fig. @J we 
calculated the same system, but we used bath with pa- 
rameters A = 5 cm" 1 , 7 = 1 ps" 1 , fiBath = 19 ps" 1 . The 
frequency is now resonant with the frequency of the elec- 
tronic coherence, which makes the effect more significant. 
The two theories give the same result for r = fs, and 
at t = 360 fs, which is approximately one period of the 
bath vibration mode. Unlike in Fig[31 the initial phase of 
pi, 12 and pn,\2 differs significantly in T = fs because 
of their different time evolution in r. 

In both the cases studied above, the time evolution of 
the off-diagonal elements of the second order response 
operator is slightly modulated by the time evolution of 
the vibrational DOF. The phase of the oscillations seems 
to be mostly unaffected. 

Figs. [S] and [5] demonstrate the influence of the reso- 
nance coupling on the population and electronic coher- 
ence dynamics in the homodimer. As above, we per- 
form calculation according to standard CL-QME and 
the parametric CL-QME. This time, we choose the over- 
damped Brownian oscillator with fixed r = 60 fs and 
A = 100 cm" 1 , reorganization energy A = 120 cm" 1 and 
correlation time r c = A -1 = 50 fs to represent the bath, 
and we change the resonance coupling. We calculate 
both diagonal and off-diagonal elements ("populations" 
and "coherences") of the operators pjiji- For J = cm" 1 , 
there is no population dynamics. By increasing the cou- 
pling, the difference between the theories in the diagonal 
elements increases. In the Fig. [5J the dipole moments of 
the molecules are anti-parallel, while in Fig. [5] they are 
parallel. 

Let us now investigate a molecular dimer coupled to 
overdamped harmonic bath and to a single harmonic 
mode with frequency fieath- The harmonic mode is as- 
sumed to continue oscillating even long after thermaliza- 
tion in the overdamped part of the bath has taken place. 
Therefore, the r-dependent term of Eq. (|5TT) changes the 
system dynamics also at long times. The time depen- 
dence of the density operator elements for different r is 
shown in Fig. [7] Parameters of the overdamped bath are 
A -1 = 100 fs and A ove rdamped = 12 cm" 1 and of the har- 
monic mode Aharmonic = 30 cm" 1 , fiBath = 60 cm" 1 . The 
dimer is characterized by A = 100 cm" 1 , J = 33.1 cm" 1 
and the parallel electronic transition dipole moments of 
the molecules. We can notice that the interaction of the 
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Figure 5: The time evolution of the real parts of 7?/ and 
Rn and their dependency on the resonance coupling j for 
a heterodimer with energy gap A = 100 cm -1 and an anti- 
parallel arrangement of the transition dipole moments, inter- 
acting with a bath represented by an overdamped Brownian 
oscillator with reorganization energy A = 120 cm -1 and cor- 
relation time t c = A -1 = 50 fs. The delay between the two 
pulses is fixed to r = 60 fs. Black lines, the matrix elements 
11, 22 (left column) and 12 (right column) of Ru according to 
standard CL-QME (dashed) and the corresponding paramet- 
ric CL-QME (full). Red lines, the matrix element 12 (right 
column) of Ri according to standard CL-QME (dashed) and 
the corresponding parametric CL-QME (full). 



vibrational mode with the electronic DOF induces oscil- 
lations in both the diagonal and off-diagonal elements of 
the density operator. The amplitude of the oscillations 
increases with increasing r. Since the EGCF of the har- 
monic mode, Eq. (1171) . is periodic, we expect the relative 
amplitude of the oscillations to decrease again for suf- 
ficiently long r, and to become zero for r = 27r/f2Bath- 
However, in such long r the second order response would 
be almost zero due to the decay of the optical coherence 
in the first interval. Fig. [Jj suggests that the phase of 
the oscillation does not change linearly with r. To un- 
derstand this behavior, we can study a simple equation 
of the form 
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Figure 6: The same as Fig. [5] but for a heterodimer with 
parallel transition dipole moment arrangement. 



p(t) = -a(p(t) - p ) + (f(t + t) - f(t))p(t), (52) 

which describes exponential decay to a limiting value po 
with the rate constant a and a modulation by the func- 
tion 



/(*) 



Z 

J dr A cos(ojt). 



(53) 



Eq. (|53")l represents a first integral of the EGCF of an 
underdamped vibrational mode, Eq. (ITTl) . The solution 
of Eq. ((5H1) witu tne parameters a = 0.01, A = 10~ 4 , 
uj^ 1 = 60 and po = 0.6 is shown in Fig. [5J We can 
see that the oscillation phase and period of the oscilla- 
tions is indeed not proportional to r, and the amplitude 
increases with r similarly to Fig. [71 The behavior is 
therefore a direct consequence of the parametric term in 
the parametric QME. 

The overall picture arising from the numerical simula- 
tions is the following: Except for special cases, such as the 
molecular homodimer, the excited state dynamics of an 
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open quantum system, as it is observed by the non-linear 
spectroscopy, indeed depends on the delay r between the 
FWM scheme. The 2D spectroscopy therefore observes a 
certain averaged dynamics. The effects seem to be rather 
small in most cases, but they might be observable by an 
advanced implementation of 2D spectroscopy. They are 
especially pronounced in the case of the intramolecular 
vibrational modes, which have frequency similar to the 
electronic energy gap between excitonic levels. Both the 
dynamics of electronic level populations and electronic 
coherences are affected. In order to identify these effects 
in the experimental data, the theory has to be extended 
to include also the third interval of the third order non- 
linear response. The corresponding projection operator 
Vt+t which now depends on the duration of both the 
coherence and the population intervals r and T, respec- 
tively, has been already proposed an tested for j mn = 
in Ref . [4j| . The formulation of the theory for the third 
interval of the response of a multilevel excitonic system 
in a similar manner as performed in this paper for the 
second interval will be the subject of our future work. 
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Figure 7: The time evolution of the real parts of Ri and Rn 
a heterodimer with energy gap A = 100 cm -1 , resonance cou- 
pling j = 33.1 cm -1 and a parallel arrangement of the transi- 
tion dipole moments, interacting with a bath represented by 
one overdamped Brownian oscillator with parameters A -1 = 
100 fs and A ovcr dam P cd = 12 cm -1 and an underdamped vi- 
bration with reorganization energy Aharmonic = 30 cm -1 and 
frequency f^Bath = 60 cm -1 . Black lines, the matrix elements 
11, 22 (left column) and 12 (right column) of Rn according to 
standard CL-QME (dashed) and the corresponding paramet- 
ric CL-QME (full). Red lines, the matrix element 12 (right 
column) of Ri according to standard CL-QME (dashed) and 
the corresponding parametric CL-QME (full). 
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Figure 8: The time evolution of "population" according to 
model Eq. (|52[) . The curves are solutions for r increasing 
from fs to 100 fs with step of 20 fs. 



VII. CONCLUSIONS 

In this paper, we have demonstrated that the third 
and the second order non-linear response of a multilevel 
system cannot be completely evaluated by propagating 
reduced density matrix by equations of motion derived 
using a single projection operator. Such treatment would 
inevitably neglect correlations between the time evolu- 
tion of the bath during the neighboring intervals of the 
non-linear response. We have derived equation of mo- 
tion, the parametric quantum master equation, which 
takes these correlations into account approximately, and 
and showed that in the absence of resonance coupling the 
method yields an agreement with the result obtained by 
the second order cumulant method. We confirm by nu- 
merical simulations that for different delays r between 
the excitation pulses, distinct dynamics of both excite 
state populations and electronic coherence occurs, in the 
presence of environmental degrees of freedom with finite 
bath correlation time and in the presence of intramolecu- 
lar vibrations. The two-dimensional Fourier transformed 
spectroscopy sees in these cases some averaged dynamics. 
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Appendix A: Third-Order Non-linear Response 
Functions 



The four standard response functions (so-called Liou- 
ville pathways) of a two-band electronic system read in 
our block formalism as 
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yields 



R x (t, T, t) = ti{p^U (e9e9) (t) V (esee) 

x l((eeee) (rp^y(eeeg)y(egeg) ^y(eggg)yy(99)} t (_^1) 

R 3 (t,T, T ) = tr{^ ge W e9e9) {t)V {e " 9) 

x U {gaaa)(T) V (sgge) U (gege)^ v {gegg) w {gg )^ ^3) 

R 4 (t,T, T ) = ti{ii {9 ^U [e9e9) {t)V [e " 9) 

x ^(99S9)(y)v(99e 9 )^(e 9e3 )^ T ^y( e3 99) w/ jg 9 )|_ 

Here, wiq 9 '' = w oq |<7)(g where w oq is the bath equilib- 
rium density operator. 

Appendix B: Derivation of Relaxation Terms 

In this section we will evaluate the term NZ 3 of Eq. 
. Applying the definitions of its component operators 
from Section W\ we obtain 



NZ, 



-tr B < 



where 

p {I) {t)=tr B [v T W { - I \t)}. 
Using Eq. (glj) yields 



(Bl) 
(B2) 



ink 



U 9 (r)UlHr)K k }w^ 
- £ e 9 s (r) tr B {p«(t) w cq U 9 {r)uf{r) 



in k 



x ArA$, 



,p m (t)}ii)W ■ (B3) 



Now we have to set e s fc ( - T ' ) » 1 since its contribution is of 
higher order in A$. Applying the first order expansion 



U°(t)U?(t) « 1 + - / rfr' A$ Sfl (-r'), (B4) 



NZ, = 



rnk v 



'cq 



x K- k ^ m {t)K m {t)p w {t) w cq jw ( /K (B5) 
Now, we define a line shape function 

t r 

ff»B(*) = J?j dT j T ' ( A$ a(r')A^), (B6) 



where it is noteworthy that one of its indices represents 
the site- and the other the exciton-basis. The NZ 3 term 
can then be expressed as 



NZ 3 = - % I dr> 



o 

X 



J2 K m {t)pV\t)K- k - pV\t)K~ k K m {t) 



jnk 



tTfl { A$fe(-T')A* m (i) Woq} 

(K m (t)pW(t)K k - pV\i)K- k K m {t) 



(B7) 



In order to obtain the last two lines of Eq. (|50|) we will 
trace Eq. ([B7j over the bath DOF. 



Appendix C: The r-dependent Term for a 
Homodimer 



Let us assume a molecular homodimer with the Hamil- 
tonian 



3 
3 



Hs-B 

This yields in the exciton basis 



A$i 
A$ 2 



-3 
j 



1 



S-B 



A$i + A$ 2 A$ 2 - A$i 
A$ 2 - A$i A$i + A$ 2 



(CI) 
(C2) 

(C3) 
(C4) 
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The derivative of the line shape function which enters 
the T-dependent correction T>(t; r) of the CL-QME then 
reads as 

* 2 

9 mk (t) = \j & J2(A$ m (t')A$ n ). (C5) 



We neglect the cross-terms (A$,(t)A$j) for i ^ j , be- 
cause we assume only local correlations of the bath. In 
a homodimer A$i = A$2, and we get g lk (t) = <7 2 fcM = 
g(t) from which it follows that 



9U(t + T)-9tt(t)=g*{t + T)-g*(t). 



(C6) 



By substituting this equation into Eq. (|51l) . the 
r— containing term T>(t\ r) turns to zero, because 



J2 {K mP {t)K- k - p(t)K- k K m ) = 0. 



(C7) 



mk 



Thus the T-dependent correction to the CL-QME is zero 
for the case of a homodimer. 



Appendix D: Pure Dephasing of an Electronic 
Coherence 



In a pure dephasing model of the system-bath inter- 
action, analytical solution for the dynamics of electronic 
coherences can be obtained with the cumulant expan- 

(2) 

sion technique. The time evolution of coherence p ee , (t; r) 



reads as 

P { 3(t,r) = tr B {u e (t)U g (r)W cq Ul,(r)Ul,(t)} pfj(0), 

(Dl) 

where we set all transition dipole moments to one. By 
using the second cumulant procedure [14j we can evaluate 
the expression into 



P^(t; T ) = e -9ee(t)-S:' e '( f + T )+9='e(t) 

e +C'(*+-)-O(-) p ( 2 J(0). (D2) 



Equivalently, we can rewrite Eq. (|D2[) in form of a master 
equation 



dt 



- g* e , e ,(t) - g ee (t) 

+ 9e'e(t)+g* e ,(t) 

-C& e '(* + r) 

+ (9L'(t + r)- 9 ;At))]p { eJ(f,r). (D3) 



The first four terms on the right hand side of Eq. (|D3[) 
are described by the CL-QME. The remaining terms can 
be matched with the term T>(t; r) derived in this paper. 
Eq. (|50p therefore gives the correct result for this simple 
exactly solvable case. 
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